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1 . Introduction 


The first part of the present Report supersedes Ref. 1; 
the second part contains a number of calculations and the discus- 
sion of their results. 

In 1976, when the second part of the work reported in [1] 
was ready to be published, I decided to reconsider the numerical 
technique exposed in [1] in order to simplify it and to increase 
its accuracy. The search was successful, thanks to the good work 
of de Neef on the methods to compute shock and body points, and 
of Zannetti on the integration scheme for interior points. The 
great advantage of de Neef’s way of computing shock and body 
points over the method outlined in [1] stems from his recasting 
of the characteristic equations in such a way that second deriva- 
tives are no longer necessary. In the present category of prob- 
lems, where the evaluation of second derivatives implies an ela- 
borate manipulation of terms connected with complicated conformal 
mappings, their elimination provides a substantial saving in cod- 
ing complexity and running time. On the other hand, the integra- 
tion scheme which took form as a result of my discussions with 
Zannetti is, to this date, the closest of all available schemes 
to the physics of the flow and it has second order accuracy; it 
is, therefore, best suited to handle complicated flows with for- 
mation of imbedded shocks, and entropy layers, in this way allow- 
ing the number of grid points to be kept at a minimum. Details 
on the shock-and-body points method can be found in [2]; details 
on the integration scheme (which I call the X - scheme) are 
given in [33. 

In order to avoid confusion, I prefer to expose the 
analysis in its entirety, rather than referring to [1] and indi- 
cating what has to be changed or eliminated. Therefore, parts of 
the first few Sections of this paper are similar, but not identi- 
cal, to corresponding parts of [1], which should be considered 
obsolete. 

As stated in the Introduction to [1], wa will provide a 
detailed description of a computational program for the evalua- 
tion of three-dimensional, supersonic, inviscid, steady flow past 
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airplanes. No imbedded shocks are considered here explicitly. 
The emphasis is put instead on how a powerful, automatic mapping 
technique is coupled to the fluid mechanical analysis in order to 
assure a high degree of accuracy without increasing the number of 
canputational nodes beyond reasonable limits. 

Care has been taken to describe and to code each of the 
three constituents of the 'analysis (body geometry, mapping tech- 
nique, and gas dynamical effects) separately, to facilitate ap- 
plications to different geometries or substitution of the present 
set of equations of motion by other sets. Sections 5 through 10 
contain the outline of the code dealing with gas dynamical ef- 
fects; all their statements and formulae are unaffected by 
changes in the mapping technique or mapping parameters or in the 
geometry of the airplane. All expressions related to the mapping 
are given in Sections 11 and 12. Sections 13 through 15 deal 
with the choice and treatment of initial conditions. Results of 
computations based on sample geometries, and discussions are con- 
tained in the remaining Sections. 


£. Frames of reference 

The free stream is assumed to be uniform, with a given 

Mach number, M^. A Cartesian, orthogonal frame of reference, 

(x,y,t) is defined as having the y and t-axes in the symmetry 

plane of the vehicle, the t-axis lying along the fuselage. The 

unit vectors of the x,y, and t-axes are called 1, 3 and K, 

respectively. The free stream velocity vector, ^ , is parallel 

00 

to the (y,t)-plane; the angle of attack, a, is the angle between 
^ and therefore, 

CO 

^ = V (3 sina + K cosa) 

00 oo 

In each cross-sectional plane, a complex variable, z, is 
as 

z = x+iy (2) 

A conformal mapping (details of which will be found in Section 


( 1 ) 

defined 
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12) defines a one-to-one correspondence between the portions of 
interest of the right-hand side of the z-plane and a portion of 
the right-hand side of a 5-plane where, by and large, the image 
of the cross-section of the airplane is nearly circular; it is 
convenient, thus, to express the complex variable 5 in the 
form: 

C = pe^® (3) 

The analytic function c(z) implies that p and 0 are functions 
of X and y, and vice versa. Such functions, in general, change 
from one cross-section to another; therefore, we may write; 

p = p(x,y,t) X = x(p,0,t) 

0 = 0{x,y,t) y = y(p,0,x) (4) 

T = t t = T 

We must take good care of denoting t by another symbol, t , when 
considered in connection with p and 3 since when t changes and 
x,y remain unchanged, p and 0 generally change; consequently, 

derivatives with respect to t (at constant x and y) generally 
differ from derivatives with respect to t (at constant p and 0 ). 
Let p= b(0,T) and p = c(8,t) be the equations of the image of the 
airplane body contour and of the image of the bow shock in the 
5-plane. A non-conformal mapping, defined by a suitable function 
of p , 0 and T ; 

X = X(p,0,t) P = p(X,Y,T) 

Y = 0 0 = Y (5) 

T = T T = T 

will transform the region of interest in the right-hand side of 
the c-plane bounded by p = b and p = c onto a rectangle, bounded 
by the lines: 

X = 0, corresponding to p = b (body) 
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1. 

corresponding 

to 

p = c (bow shock) 

(6) 

-ir/2 , 

corresponding 

to 

0 = — ir/2 (windward 

symmetry line) 

ir/2. 

corresponding 

to 

0 = ir/2 (leeward 

symmetry line) 


An example of such a function, X(p,e,t) will be discussed in Sec- 
tion 11. 


Derivatives related to the mappings . 


Let 


g = 

® dz 


Ge 


lu 


(7) 


be the complex derivative of c with respect to z (at t, t ,T 
all constant); similarly, let 


- 1 d log g 
“ g dz 


4'^+i<t'2 


From (3) and (7) it follows that 


where 


51 - e/g _ gi(0-m) 
pg “ mod(t/g) 


l2+iS 


( 8 ) 


(9) 


t - cos(6-(i)) , S = sin(e-to) 


( 10 ) 


We introduce now the notations, ip and f, for two analytic func- 
tions obtained by differentiating g and c with respect to t 
(that is, at constant x and y); 




f = a l og c = f + , f 
^ 3t 1 ^ ^ ^2 


( 11 ) 

( 12 ) 
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Recalling that 



Derivatives related to the mappings 


de _ a( pcose) . a( psine) _ a( pcoso) 3( psin9) 
dz “ 3x ^ ^ 3x " 1 3y ^ 3y 

and 

31og e _ 1 M • 39 
3t “ p 3t ^ ^ 3t 

we obtain: 


Pjj = GU. Py = G2, 

6 = -- 2, 9 = -e, e^ = f 
X p ’ y p ’ t 2 


X = 0, X = 0, T. = 1 

X ’ y ’ t 


Conversely, noting that 




we obtain: 


X = C/G, X- = - pS/G, X = - (jCf,-Sf_)p/G 
P u T I ^ 

y = S/G, yg = pC/G, = - (2f^+2f2)P/G 


t = 0 . t = 0 , 
P 9 


t = 1 

X 


Between the two sets, (p,6,x) and (X,Y,T), the following 
tions hold: 


Px = 1/x^. Px = - X/X^. Pj = - x^/x^ 

9y = 0, 9y = 1, = 0 


Xy = 0, Ty = 0, Xy = 1 


^p “ ^0 ■ " ‘’y'^^X, V " ” ^T^^^X 

Y = 0 , Y. = 1 . Y = 0 

p 9 X 


Tp = 0 , To = 0 , = 1 


By combining (17) and (18), we obtain: 


( 13 ) 


(1^) 


(15) 


(16) 


(17) 


rela- 


(18) 


( 19 ) 
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P /- 0 




' GX ’ " "G ^^pX ■ " G X ■ ^^^1 “ ^^2) G 

P P P 


^X ■ GX ’ " G pX ~ “ G X “ ^^^1 ^^2^G 


t 


tx-0 


4 = 0 


. t„ = 1 


( 20 ) 


The following formulae are also obtained easiiy: 


Gp = G(j)^/p , Gg = - G(t>2 , G^ = G Ci|i^-<t)^f ^+4>2f 2^ 

Up = <>2^P » “0 “ “**1 * “t ” 

G^ = 4 ~ 0^(<5(J)^-l2<t)2^/p» 0^^ = G<;)^ 

= “ G(S())^-l2(^2^/p* Uy = G(£(ji^+2(t)2)/p, oi^ = ^ 
4. Important unit vectors 


(21) 


( 22 ) 


We begin this Section by defining a p-line on a 
t=constant (physical) cross-sectional plane as a line along which 
0= constant; similarly, a e-line will be a line on the t = con- 
stant plane along which p= constant. The unit vectors, i and j 
will be used to identify the tangents to a p-line and to a 0-line 
respectively. Note that 


i = ei + 2J , 1 = 2i - S3 

3 = - !Si + SJ , 3 = S3 + U3 

By using (20) and (23) for any point, 

Q = xl + y3 + tR 


(23) 


(24) 


we obtain: 
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Q 


Y 



Xa 


J) 


p 

■ G 


* f,) 


1 - 


G y 


( 25 ) 


The unit vector, R, normal to an X = constant surface, is 
important for the calculation of body and bow shock points. The 
body, indeed, is defined by X = 0 and the bow shock by X = 1. In 
general. 


f) 


N^i + NgJ + = 


mod(Qy><Q^) 


(26) 


where Q is a point on the surface. From (25) it follows that 


1 

= V • ^2 "" PX- ’ "3 


= N^d 


(27) 


with 




'Xx Xg ' 

Xp XpV‘’^1. 


V - 


Xe 2 2 

1 . 


1/2 


J 


(28) 


In particular, at the body, from (18), 


Xg/PXp = - b^/b, X^/Xp = - b^ (body) (29) 

Note also that we can write b^ or bg, and b^ or b^, indifferent- 
ly. Therefore, at the body (27) and (28) take on the form; 


d = 


N^ = 1/'’, 

(b^+bYf2-bf^)/G, 


= - (bY/b)N^ = N^d (body) ( 30 ) 

= Cl + (b^/b)^+ (body) ( 31 ) 
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Similarly, at the shock, 

X„/pX^ = - Cy/c, X^/X = - C-. (shock) (32) 

op I TP 1 

= l/v, ^2 ” “ = N^d (shock) (33) 

d = - (c^+CYfg-cf-j)/^, V = [1 + {.Q^/c)^+ (shock) (34) 

Let 


F(x,y,t) = 0 


(35) 


define the geometry of the body in the physical space. The image 
of the body in the (p,0,t) space is 


p = b(0,T) 


(36) 


To evaluate (30) and (31), that is the normal to the body, we 

need b /b and b . At x= constant, 

0 T 


* * 6 > * ^ ° 


(37) 


Consequently, and using (17): 


b “ 2F + SF 
X y 


(38) 


Similarly, at 0=constant, 


F(xb +x)+F(yb +y)+F. =0 
xpx X y pT X t 


(39) 


and 


^x " " ^2*^0 + Q UP + 2F 


(40) 


Therefore, at the body. 
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d 


F. 



+ SF 

y 


(body) (41) 


and (30), (3D can be replaced by the simpler expressions: 


= (eF^+2Fy)/v, Ng = (CFy - ZF^)/v, N, = F,/v 


V = (F^ ^ . F^) 


2 , 1/2 


’3 " * t' 
(body) 


(42) 


The above formulae are general. For any particular geometry, 

F , F and F must be evaluated. 

X y t 


Equations of motion 


Having chosen a suitable reference length, the 
pressure, density and temperature of the free stream are chosen 
as reference pressure, density and temperature, respectively 


(P. 


ref 


"ref =■”“ «ref> ’ 


quantities, the equation of state is then 


With p, p, 0 measuring non-dimensional 


p = pe (43) 

The reference velocity, is defined by 

“ref “ Pref'^^ref “ ^®ref 

where R is the gas constant divided by the molecular weight of 
air. The speed of sound in the free stream, in a non-dimensional 
form, is then 


= /y (45) 

The logarithm of pressure is denoted by P: 

P = In p (46) 

A non-dimensional entropy, S (which is the difference between the 
local entropy and the free-str.eam entropy divided by c^) is re- 
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lated to non-dimensional temperature and pressure by 

S r Y In 0 - (y- 1)P , e = exp[(Y-l)P/Y + S/y] (47) 
Euler's equations of motion in non-dimensional form are: 

V.7P + yV.V = 0 

- ^xVx^ + 0VP = 0 (48) 

^.VS = 0 

With R = K, let 

i? = w(x + K) (49) 


where 


X = oi + nJ 


(50) 


and let 


V, = f- I + I- ^ 

1 3x 3y 


(51) 


Note that 


^.7P = w(x+K) .(V^P+Pj.{<) = wx.V^P + wP^. 

/>- 

V.V = , [w(x+J<)]+Mj^ = X.7 ^W+w7^ .X+W^ 

„ /-2 1 2 „ -*-2 1 2 .-^ 2 . ^ 

— V(V ) = wVw(x +D+ ^ 7^X + ^ (x )|.K 

^xVx? = - w(7^W.X+M^) (X+^f)+w7w(x +1 )+w'^[(x^.x){c-X^.+Xx7^xx] 

therefore, (48) takes the form: 


(52) 


w(x.V^P+P^)+y(X.V^W+wV^ . x+w^) = 0 
x^+w( w . X+Wj^) X+W^Xj^-W^Xx xx+0 P 


w(7^w.x+w^.)+0P^.=O 


= 0 


(53) 
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X.7,S + = 0 

I O 

The third of these equations can be used to simplify the first 
and second equation; finally, the following system is obtained: 


(1 - + x.v.p + yv .X = 0 

w 


1 -►2 -► 

- X’^ViXx + — (V,P - Pt_x) + Xv = 0 


2'’r 


w 


(54) 


X.V^S + = 0 

The third of (53) is not needed; the above system is composed of 
four scalar equations for the two unknown scalars, P and S and 
the two-component unknown vector, x- Once P and S are deter- 
mined, 0 is obtained from (47); the modulus of the velocity, q, 
is obtained from 


- ^X_ (0 _ 0) (55) 

Y“1 o 

where 0^ is the (non-dimensional) stagnation temperature, and w 
follows from 


2/1 2 2 . 2 faei\ 

w(1+a +n)=q (56) 

There are definite advantages in using (54) as a basic system of 
equations (instead of (48) or of equations in divergence form). 
It contains only four differential equations to be integrated, 
and it provides a clear separation of unknowns, S on one side and 
P and X on the other side, which is particularly welcome in prob- 
lems where strong entropy gradients occur X4,5]. Another advan- 
tage of (54) stems from the fact that operates on the 

(x,y)-plane only; therefore, it can be expressed in terms of p 
and 0 as independent variables, using i and j as unit vectors. 
In particular, note that 

'l'’ = ^ } ''o3) 
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f,2 

X^V^xx = — C(;i^) - (^) ](ni-oj) (57) 

I pup U 0 

’r* = 

= GC(oo^+nnp)i+^(crOg+nng)j3 
x^ = ■'■ 

These expressions can be substituted into (54); in doing it, how- 
ever, note that t also must be substituted by t, and that, for 
any function, 

$ - $ + ® p. + $-9. = $ + pf,4 + f^$. (58) 

t T p^t 9t X Ip 20 

Using the notations: 

, a2 

w 

A, = G- + pf, , A„ = ^ + f^ 

1 1C 1 2 pK 2 

(59) 

= Go + pf^ , ■ ®2 " ^ ^2 

p 

D = — Cn(1-<ji-.) - Oijio] ■*■ ^2 “ '**2 

^ (J * £- fc, fc. 

and taking (21) and (22) into account, (54) become; 

P +A.P +A„P + ^(o + ^ )+ Ca(1-4>. )+n<i>o3 = G 

xlp20< pp0 icp 1 2 

0 + 6 ^ 0 + 8200 + ^[-oP^+(G-opf^)P - ofgPeJ-nD = 0 
w 

( 60 ) 

n +B,n +B„n„+ -^C-nP "PtfiP„+(7 -nfp)P„ 3 +oD = 0 

TlptO d T Ipp.^U 

W 

^x'^®l^p‘*'®2^0 " ° 
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The final form of the equations of motion is obtained by 
expressing the derivatives in terms of X, Y and T, considering 
that, for any function 4: 


4 = 4^ X 

P X p 


^ " *Y *x 


(61) 


4 = 4- + 4y X 

T r X T 


and, consequently: 


where 


4 + B.4 + B„4„ = 4™. + E4v + B„4v 

X Ip 2e T X 2Y 


4 + A.4 + A„4„ = 4™ + C4y + A„4v 

T Ip 20 1 X 2Y 


(62) 


C = X +A,X +A_X„ , 

T 1 p 2 0 


E = X +B,X +B-X„ 

X 1 p 2 0 


(63) 


With the additional notations: 


L = a(1 — 4>-j ) ■*■ 42 


F = - aX + (G-apf,)X - af_X„ 
X Ip 2 0 


(64) 


H = - nX^ - pnf^X^ + (7 - nfp)X. 

X Ip p ti a 

the equations to be integrated at every grid point are: 

f 1 (n,*x,nx U] = 0 

a.p+Eox+B2a^+ ^ [-aP^+FP^-af^P.^] - nD = 0 


’^T^^^X^^2^Y'*' ^ [-nf’.j«+HPjj+(^ -nf’2^^Y^ + oD = 0 
w 


( 65 ) 


S^+ESj+B^Sy = 0 
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6. The integration method 


The equations of motion are integrated using a 

predictor-corrector method. If we denote by f the initial value 

of P, 0 , n or S, by ? the corresponding value at the end of the 

N 

predictor stage, and by f the value at the end of the corrector 
stage, that is, the final value, we use the formulae: 

? = f + f^AT (66) 

and 

f^ = |-(f + ? + ?^AT) or f^ = ?+j(T^-f^)AT (67) 

in the predictor and the corrector stage, respectively. 

Many different integration schemes can be devised, all 
making use of (66) and (67), the difference residing in the way 
f^ and are defined. For example, the original scheme suggest- 
ed by HacCormack [6] used the equations of motion in divergence 
form and discretized the space-like derivatives using forward 
differences in the predictor and backward differences in the 
corrector (or vice versa) . Without recasting the equations in 
divergence form, I used the same alternating forward-backward 
differencing in a large number of works, and I called it the Hac- 
Cormack scheme. The entropy equation, however, was always in- 
tegrated approximating the derivatives with upwind differences, 
in order to maintain consistency of the numerical approximation 
with the physical nature of the problem [4]. 

Numerical work on three-dimensional, steady, supersonic 
flows which began under promising auspices, showed that, as the 
body geometries became more complicated' and the local Mach 
numbers closer to 1 , the MacCormack scheme was losing accuracy 
and reliability. Since the most probable cause of failure was a 
lack of consistency between the physical domain of dependence of 
a point and its numerical domain of dependence, in the second 
phase of the work related to the present paper I tried to approx- 


- 14 - 



The integration method 


imate Lagrangian derivatives oy upwind approximations, and I at- 
tempted to give an overall second order accuracy to the scheme by 
a proper definition of difference operators in the predictor and 
corrector stages. Finally, I decided to adopt the \-scheme C3J 
because it seems to offer a maximum of domain-of-dependence con- 
sistency together with second-order accuracy. Referring to [33 
for a discussion and a detailed description of the scheme, we can 
see how the scheme is applied to the present problems, in what 
follows. 


terns; 


The first three equations (65) are split into two sys- 


and 


<=‘“x*V 


yGX, 


rGL 


+ Tiv + 

A KP X icp 


=0 


X ^ ^ o 

®T 2 
w 


Ea, 


- tiD - Pjj = 0 
w 


( 68 ) 


* M P 

hf + 2 ^x 

w 


+ Etiy + OD - ^ = 0 


P^ + 


+ T1 

^ KP ’’x 


fyPy + Bp 
^ w"^ Y 


6a pY 

2 ^T 
w 


= 0 


= 0 


(69) 


y 9 /G « _ 0n dX ri 

Hf + ~2o * ®2'’Y 2 ^T ” ^ 

vt‘^ ^ w'^ 


Using the notations of (33, 


a -c a -X5.x a s^Gv c=x5k 

®11 ■ ^ , a^2 - K ^p* ^13 - KP ^9’ KP 


i ^ a - E 
21 " 2 ’ 22 ~ 
w 


0a 


ap.. =0 , c = -tiD, kp = ^ 

^ (70) 


0H „ 

^31 ^ ^2’ ^32 ■ ° 


^33 = ^22 = °3 = ^ 
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b^1 = A 2 


bi2 = 0 


b = ^ 
13 KP 


21 


9g 

w2 


^ 2 ’ *^22 " ® 2 ’ ‘^ 23 "° ’ ^2 ~ ^2 

w 


( 71 ) 


^31 ' 2 ^p “ ’’^ 2 ^’ ‘’ 32 "°’ *^33 ■ *^22 " ® 2 ’ *^3 " ^ 

w w 


Therefore, 


‘11 


+ a 


22 


+ + ^ 12^2 ~ 


and, consequently. 


= C(aX +nX /p)^ + <(X^+X^/p^)]^^^ 

pa P 0 


= (n^+ic)^^^(72) 


xX _ c - ^ 

^1 ■ ^ KW ® 

» aG Y 
^1 "" ^2 " Kwp ® 


4 = . 

A aG „Y ' 
X„ = A_ + B 

2 2 KWp 


X^ = E 
^3 " ^2 


= a3,-xf.a,2V^3‘'3 = 2C-a„-X* = C-x^ = |2 

= »22-»r'x,3‘'3 = 2A2-b„-x![ = Aj-X^ = ^ b'^ 


(73) 

(74) 

(75) 

(76) 


X Y 

To determine and P^, Eq. (51) from [ 3 ] and its Y-counterpart 
are used. Similarly, to determine the auxiliary quantities 
it>^. '1'^ and > 1 (^, Eqs. (52) and (54) from [ 3 ] and their Y- 
counterparts are used. Moreover, 



®12 

"13 



0 

*^13 

t> 

X 

II 

X 

“32 

X 

^33 

t 

a" = 

Y 

^32 

Y 

^33 



®13 


- 


/ b 
9 °13 

0? = 1/a" 

/ 

X 

^33 

• 4 

= 1/A^ 


,Y Y 

^ ^33 


(77) 


( 76 ) 
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a J 

ai2 


0 

Tip = 1 /A 


Y Y 



X ,x 


,Y 


U32 ♦ 


U32 


Once expressions of the type and have been evaluated, is 
obtained as 

(80) 

All X- and Y-derivatives are approximated by expressions, 
generally denoted by fYi» ^'y 2’ s'”® defined in 
Section 3 of [ 3 ]. 

Approximations to and are always defined using 
upwind information only; two- and three-point formulae of the 
type shown in (14) and (15) of [ 3 ] are used in order to provide 
second-order accuracy. 


7 . Characteristic equation for body and shock points 


Although the following discussion could be conducted on 
the basis of the characteristic equations of the preceding Sec- 
tion, it is simpler to proceed without splitting the equations, 
as follows. We build up a characteristic equation in the (X,T) 
plane by multiplying the first three equations (65) by 
, U 2 , and y^, respectively, adding, and calling X the slope of 
the characteristic in the (X,T) plane: 

0 

Cy^^ — ^(ay2+nyo)3(Pj+XP^)+y2(o^+lOj^)+y2(hp+Xn^) = 

w ‘ ( 81 ) 

where R contains all Y-derivatives and all non-differentiated 
0 

terms. In turn, X is defined by 
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C-\ -^(F+aX) 

\r 


•^(H+nX) 

\r 


that is. 




E-X 


0 


0 


= 0 


E-X 


( 82 ) 


(E-XKC-X) - ^[X (H+nX)/p + X (F+oX)] = 0 
.^2^ 0 P 


. or 


2 

X^ - 2(X +A,X +A„X^)X + EC - ^(FX +HX„/p) = 0 
T 1 p 2 9 ^2^ p 0 


which, after some manipulations, yields: 


X = C±B— 

WK 

B = C(aX +nX./p)^+K(X^+X^/p^)J^^^ 

p u p o 


(83) 

(84) 


The lower sign and the upper sign must be used at body and bow 
shock points, respectively. From (82), we obtain 


u -E— X u - - y u - — X 

y^-h X, y^- ^ Xp, y^- Xg 

0 A \ X/, 

-“I + -2<»‘'2"-™3’ = 
w 


(85) 


( 86 ) 


The compatibility equation (81) is then 


±(P^+XP^)+^[Xp(a^+Xaj^)+Xg (n^+Xn^)/p] = 


(87) 
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Treatment of body points 


At body points, the boundary condition, 

1 . n = 0 ( 88 ) 


yields 


a - nby/b + d = 0 (89) 

if (M9)f (50), (30) and (31) are taken into account. On the oth- 
er hand, (63) and (29) give 

E = X (-b„ + B, - B„b„) • (90) 

p 1 I ti i 

and, using (59) and (31), it is easy to see that 

E = 0 (91) 

Following de Neef's suggestions [2], we evaluate the 
pressure at the body as follows. If we use the formulae of the 
preceding Section to obtain P^, and at a body point, we ob- 
tain three values which satisfy (87) but do not necessarily 
satisfy the boundary conditions. Let us denote them by a super- 
script E. On the other hand, there exist a set of values, 

and n.j. which are the exact solutions and consequently 
satisfy both (87) and the boundary condition. If we write (87) 
twice, once for the E-values and once for the exact values, and 
subtract, we obtain 


?T - Pf - '’T“°T‘^b 

Considering now that, across the interval AT, it is, in general, 
f^ = (f-f^)/AT, f^ = (f^-f^)ZAT (93) 

R 

if f^ are initial values, f are exactly updated values and f are 
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values updated following the procedure of the preceding Section. 
Therefore, (92) can be substituted by 


P 



^ Is 



E E, 
a +— n ] 


(94) 


Using (89) and the first of (3D* the equation to determine P at 
a body point is finally put into the form: 


f Is ^ n^(bf,-t.jf2-b^)/G] (95) 

The kinematical unknowns at body points can be determined 
as follov/s. Let v be the velocity component tangent to the body 
in the cross-sectional plane: 

V = w (n + oby/b) , (96) 

The third of (53) in the (X,Y,T) frame, with E=0, reads: 

w„+B„w„+ ^ [P„ + (X^+pX„f,+X-f-)PY + f_Py] = 0 (97) 

T2Yw T T P102X 2Y 

~ 2 

If (97) is multiplied by v/w and added to the third of (65) and 
the second of (65) multiplied by b^/b, the following equation is 
obtained: 


Vt — B2V « ta(^)^+aB2(^)Y - ^ 

pw (98) 

Note that the Lagrangian derivative of v, as expressed by (98) in 
the (X,Y,T) frame, depends on the geometry of the body and on the 
Y-derivative of P only. It is crucial to approximate v^ using 
upwind information C5J. 

From the body geometry and local values of v and q, the 
corresponding values of a and n are obtained as follows. First, 
V is evaluated using (3D; than, from ( 96 ), (56) and (75): 

w = — [q (1 + by/b ) - V ] (99) 
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(v/w)(bY/b) - d 

1 + (by/b)2 (100) 


n = v/w - o b^/b (101) 

9_. Treatment of bow shock points 


Let 


= N^i + NgJ + (102) 

be the unit vector normal to the bow shock surface. The values 
■of are the same as the values defined by (33). The 
velocity component normal to the shock in front of it, u^, is 


u = ^ = u M, + 

OO 00 00 ^ 


V N., + w N, 


(103) 


where 


u = V 2 sina 
00 00 

V = V 0 sina (104) 

00 00 

w = V cosa 
00 00 

If"we denote by u the corresponding velocity component 
behind the shock, the velocity vector, V behind the shock is; 

^ ^ + (u - u )R (105) 

00 00 


The Rankine-Hugoniot conditions provide the increment in P and 

the ratio u/u across the shock (here P is the logarithm of pres- 
00 

sure behind the shock; let us keep in mind that P_^=0): 


p 

,2 , /~2 Y-1•^ 

= In — - + In (u - -hr-) 

Y+1 ” 2 

Y-1 - 2y 1 

U = - — r U + — *-:r — 

Y+1 ® Y+1 U 

(106) 


CO 

(107) 

Since (106) and (107) 

are identically satisfied at any T, 

and P 
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and u are functions of u alone. 


Where 


P =^u 
T 3u «T 


3U 

“t ■ “»T 


2u 


3P 


3u -2 

09 u - 


3u _ Y-1 2 y 1 

1 * 3u_ Y+1 Y+1 ~2 


In turn, from (103), (9), (18) and (21) 


and 


^ooT “ ” 


'^ooX “ ‘^oo'^X 


w„x = ° 


“t " '*'2”^2'*’l“^l‘*’2'^‘*’2°T^° 


Therefore, 


( 108 ) 

(109) 


( 110 ) 


( 111 ) 


( 112 ) 


V = -(v„N^-u„N2)o)^ + + v„N 2^ + w^N^^ (113) 

We can proceed now as at body points, by writing first a charac- 
teristic equation similar to (91 ), as obtainable from the in- 
tegration of Euler's equations as executed in Section 6; 

Q 

4 ^ S - -;r^4 + ^hy)] = (114) 

T APapT XcT.X o 

where, as in the preceding Section, the superscript E means that 
the values are compatible with Euler's equations but not neces- 
sarily with the Rankine-Hugoniot conditions. On the other hand, 
the same equations can be rewritten for the exact values, which 
satisfy the Rankine-Hugoniot conditions: 
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?! + Ba ^p*-°T “ c ■*■ " ^o 


( 115 ) 


By subtracting (114) from (115), we obtain; 


?! - Pf ■*■ “ °r ” c ~ 4j)] = 0 


(116) 


Now, P^, 0 ^ and can be interpreted as the sum of two terms, as 
follows; 

f - f* + c 
T ■ T 3c„ TT 


if is obtained from the Rankins-Hugoniot conditions under the 
assumption that c^ remains constant. From (34), 


3d/3c^ = -1/G = 


3v/3c^ = dCg/v =C^ 


(117) 

(118) 


From (33), 


3N^/3c^ = -C^/v = Cg 


SN^/ac^ = -CyCg/c = Cg 


3N3/3c^ = Cgd^N^C^ = 


(119) 

(120) 
(121) 


Then, 


3u /3c_ = u C^+v Co+w 

o> X <x>5ooQo>10 12 


( 122 ) 


From (108) and (109), 


3P/3C^ --it ^12- ^4 . 


(123) 


3u^/3c^ = C^2 = 


( 124 ) 


Then , 
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k 

With 

3(u-u„)/3c^ = C^g-C ^2 = ^18 

(125) 



V = ui + vj + wft 

(126) 


it follows 

from (105) that 




u = + (u-5 )N, 

OD CD 1 

V = + (u-Q )N- 

09 CO J 

w = w + (u-u )N_ 
00 00 < 

(127) 


• Therefore, 

3u/3Ct = + = C^o 




3v/ac^ = = 0^2 

(128) 


and 

3w/3c^ = C^gN 3 +(u-a„)C^Q = 




3a/3c,p = (C2 q“OC2jj)/w = ^26 

(129) 



3ri/3c,p — C ) / w — ^20 



By 

replacing (123) and (129) into (116), the 

latter be- 


comes 



- 

ia *-°T'^^26°Tr"''T ” c ^ ’^T'^^28°TT"’^ 

p)] = 0(130) 


Let 



- 

Then, 

“7 = - "‘Si. " b! *p‘'=26 - r <=2S” 

(131) 
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°TT ” Ba ^°|"°T ~ c (132) 

Using the same argument as in (93), the increment in c^ across 
the entire integration step, Ac^, turns out to be: 

Ac^ = D^(P»-P^+ II Xp[o»-a^ - -^(n^-n^)]) (133) 

The value of c^ can thus be updated and, by a further in- 
tegration, c itself can be updated. Once a new shock geometry, 
p=c(0,t), is obtained, the new fl and u^ are evaluated; then, 
(106) and (107) can be applied to compute P and u . The three 
velocity components follow from (127); therefore a (=u/w) and 
n =(v/w) are made known. Finally, S is given by 

S = P - Y ln(u /u) • (134) 

00 

10 . General outline of one integration step 


The equations obtained in the preceding Sections are used 
to proceed from a station, t, to a station, t+At, as follows. 

Predictor stage 

Given, original values of P, a, n, S, c, c^, c^, b, b^, 

b_, c, z,’ G, f, (b and 41 at all nodal points (including body and 

^ 2 

shock points), compute q ,0,w and all auxiliary quantitites 

necessary to determine and X^, (i=1,2,3), where X^ and X^ are 

'streamline slopes' used to identify the upwind direction for the 
proper approximation of X- and Y-derivatives, when needed; 
specifically, X2=E and X^sB^. Compute all X- and Y-derivatives 
Compute P^,a^,n.p,S^. At all body points, compute v^. At all no- 
dal points, update P,a,n and S using (66). At all body points, 
update v using (66). At all shock points, update c using the 
formula : 

0*^ = 0+ c^AT + j c^^AT^ (135) 

At all shock points, compute c^ (from the updated shock 

geometry), d ,v ,R ,u^,v^,u^,P* ,u* , all the coefficients needed to 
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evaluate D^, Ac^, and an intermediate value of c^, 

Cf = Ct + Ac^ (136) 

Move to the station defined by t+At, compute the new 
geometry of the body and b^. 

Corrector stage 

The computation is restarted as at the beginning of the 
predictor stage. The geometry and the grid, however, are now 
those of the new station, and all the variables have their 
predicted values at t+At as well. The updating of P,a,ntS and v 
is performed using (67). 

At all shock points, d is recomputed using the 
value of c^, and all the other values mentioned in the 
stage are recomputed. An equation similar to (136) is 
determine the final value of c^: 

c^ = + Ac^. (137) 

where c^ is, once more, the original value at the beginning of 
the step. Once c^ is updated, a new value of d is obtained from: 

Ad = Ac^/G (138) 

with the new value of d, v,{}, and u_^ are determined; then one 
proceeds as mentioned at the end of the preceding Section. 

At all body points, (95) is applied to update P; then, 

2 

from P and the updated value of S, 9 and q are obtained; a and n 
follow from (99), (100) and (101). 


original 
predictor 
used to 
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11. Coordinate normalization and grid stretching 


The object of the transformation (5) is twofold. It de- 
fines a variable X which is constant (equal to zero) along the 
body and also constant (equal to 1) along the bow shock. In ad- 
dition, it provides a stretching of coordinates according to 
which evenly spaced grid points on the X-axis correspond to 
unevenly spaced points on the p-lines. The latter property is 
used to accumulate 0-lines in the vicinity of the body where a 
stronger resolution is needed. 

The values of the derivatives, X , X„ and X , depend on 

’ p’ 0 t’ 

the choice of the stretching function X(p,0,t). In this Section 
we give an example of such a function. If the definition of 
X(p,0,t) is changed, the definitions of X , X^ and X^ must be 
changed accordingly. The rest of the program does not need to be 
altered. 


Let 


p = c+0 tanh[a(X-1)] 


(139) 


with 


0 = 


c-b 

tanha 


(140) 


Obviously, X as defined by (139) satisfies conditions (6). Dif- 
ferent values of a provide different degrees of accumulation of 
0-lines near the body. To give an idea of the effect of a, let 
b=0, c=1. Fig. 1 plots X vs. p for various values of a. Clear- 
ly, strong stretching effects begin to appear for values of a 
larger than 2. 


It is easily proven that; 


X = 


P aC0‘^-(p-c)'^] 


P-c 


(c^-b„)+0c^]X„/0 


^0 "*-tanha'''Y"^Y'^‘'^Y"“p' 


(141) 

(142) 
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= -r.7f^(c.-b„)+0c_]X /0 
T tanha IT T p 

(143) 

On those rare occasions when stretching 
the following definition of X can be used: 

is not needed. 

x-”-? 

c-b 

from which 

(144) 

X = -^ 
p c-b 

(145) 

X„ = [(X-1)b„-Xc„]X 

0 lip 

(146) 

X = [(X-1)b -Xc„3X 

T lip 

(147) 


12 . Explicit computation of terms related to the mappings 


The mapping of the z-plane onto the c-plane, mentioned in 
Section 2, is performed according to the general scheme exposed 
in [73. In the z-plane, the ’hinge-points' are denoted by h^^ 
(1=1 through J). In addition, hj^^ ^ and hj^g i ^re the affixes 
of the lower and upper intersection of the cross-sectional con- 
tour with the y-axis. 

Let 

z^ = z (148) 

and a sequence of J mappings be used, each defined by the equa- 
tion: 


^J+1 


-6 .o . 


■1 .1 _ 
+ 6 a ~ 
3 3 


(z.-h ) 

J JJ 


z +h * , 
J JJj 


6 

3 


or its inverse: 


( j=1 ,2,.,.,J) 


(149) 


z .-h . . 

_J JJ _ 

z +h * ” 
J JJ 


(z. a.)l/6 

J+1 J J 

z .+6 a . , 

J+1 J J; 


( 150 ) 
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Fig. 1 


where 


hj. = * 1 (151) 

The mappings are automatically performed in order of increasing 

With 



J+1 I J+1 




(152) 


the t-plane is defined as 


t = Vl - 

To obtain the derivatives used in the c(?mputation of the 
flow field, the following definitions are used: 


k 


1 
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2 2 2 
k- = zf ,-6 
2 J+1 J J 

'=3 = -“j"! 

'=4 = 

kg = log (kg/kj^) 

# 

k^ = o. 

6 J 


"7 = 

kj = Zj-iBj 
«5 = 

k,o= aj5j(aj5j+».5j) 

k^^= 1 /(kgk^) 


v;here dots mean partial differentiations with respect to 
constant x and y (that is, at constant z) . We also define 


®j ■ dz. " 
J 


and we note that 


3 z. . . 3 z . . 

■ ■ h + ■— - h.^ = g (k,k,-h. .) 

3hjj jj 3hj* jj *^3 16 jj 


Therefore , 


= ej(k^4k,k5)4[Zj^,k,/a^SjH-5k2k56^)/aj6j 


and 


• • 


• • 


^ “ ^j+1" 2^'^j+i,j+i'^'^j+2,j+i^ " ^j+r^ 


( 154 ) 


t at 


( 155 ) 


( 156 ) 


( 157 ) 


( 158 ) 
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The t-derivatives of the h^^ can be obtained as particu- 
lar cases of (157)# with h^ h^j in lieu of # Zj respec- 
tively. In particular, if i=j, h. . =6.o., as evident from 

089). and 

The functions g, <l>, f and tl) are determined at each grid 
point as follows; 


!! j±i i • 

® ■ ^^1 “ j=i 


(159) 


♦ = 


- 1 d log g _ ^ 


J d log g. ^ J d log g. j-1 


A - 1 


g dz 


8 9=1 


I [- 


^ 9=1 ■ i="o 


® j =1 2 1=0 


(160) 


with g^=1 , 


? at c '■^j+1 


( 161 ) 


^ = U^S_g = 1 b. 

3t , g. 

J=1 J 


(162) 




At every new cross-section, the geometry defines the 
values of h^.j , h^.j (1=1 through J+2) as well as of 6 ^ , 6.) (j=1 
through J). According to the procedure explained in' ’ ni the 
mappings are performed in order of increasing 6 ^ . First, the 
values of h^^ (j=1 through J) are found by repeated applica- 

tions of tlie mapping routine, which also provides the values of 
3h, . ,/3h, .. Then, 2 is obtained from (152). Repeated applica- 
tions of (157) allow all the values of h^^ to be determined. It 
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is possible, thus, to evaluate £ as well. 

Grid points on the contour of the body are determined 
next. For them, the value of 0 (that is, Y) is assigned and the 
value of p (that is, b) has to be found. Letting aside those 
simple problems (such as circular or elliptic cones and cambered 
wings as defined in this Report) where the body contour can be 
mapped on a circle with a known radius, b, the task can be per- 
formed in two different ways. The first seems to be more direct 
and, consequently, less time consuming, and it proceeds as fol- 
lows. 


Given a value of Y, a value of b is guessed. All the 
mappings are executed in reverse order until a value of z is 
found. Using its x, the geometry subroutine is called to gen- 
erate the corresponding value of y. Then, the modulus of z and 
the modulus of x+iy are compared and a trial- and-error procedure 
is started in order to get coincidence of the two moduli. 

Tne procedure, if working at all, works very fast; in 
fact, the first guess for b can be taken from the final result of 
the preceding point and it must be very close to the correct 
value, since the body contour in the c-plane is almost circular. 
Unfortunately, certain geometries require the use of tri- 
gonometric functions and square roots in order to get y, once x 
is given. In the trial-and-error procedure, values of x which 
produce, for example, negative radicands or sines greater than 1 
may appear. A typical case is shown in Fig. 2 where the abscissa 
of point z^ is able to provide a body point, z'^ but the abscissa 
of point 7 .^ cannot generate any body point. Round-off errors 
which vary from one computer to another compound the difficulty. 
Constructing a procedure which takes care of such possibilities 
and corrects the guesses in such a way that a body point is al- 
ways definable in the course of the trial-and-error iteration, 
and which also provides a final value of z with the required ac- 
curacy seems to be a rather difficult task. • For this reason, the 
procedure described above was abandoned in favor of the following 
one, whose fine details seemed to be easier to implement. 

Given a value of Y, a value of x is guessed (suitably lo- 
cated between the value of the previous point and the maximum 
value admissible at a given station) . From the geometry 
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subroutine, a corresponding value of y can thus be found, with no 
exceptions. By applying the mapping procedure to z (=x+iy) , a 
value of 5 is found, and its argument is compared with the given 
Y. the procedure is repeated by increasing or decreasing x until 
two successive values of the argument of 5 are found which brack- 
et the given Y; then, a trial-and-error iteration is performed. 

Subsequently, g, 41, 41, f, E and S are computed and the 
geometry subroutine is used to determine F^, and F^ at all 
body points. 

Once b has been determined and c is kno;m from the updat- 
ing of the shock, all grid points in the t-plane can be defined 
according to ( 139 ) and ( 140 ). For each of them, the mapping pro- 
cedure is applied in reverse order until each grid point in the 
physical plane is obtained. In so doing, we also obtain all 
values necessary to evaluate g, fi iC and S. 


13. Initial conditions 


Three different assumptions can be made with respect to 
the form of the front end of the vehicle: 


- 33 - 



Initial conditions 


1) The body begins with a blunted nose, or 

2) the body begins with a pointed cone, or 

3) The body begins with a hollow intake, having a sharp 
lip contained in a plane normal to the t-axis. 

The first assumption is not only realistic since most of 
the vehicles have a blunt nose, but is also convenient for a safe 
start of a numerical analysis since the three-dimensional blunt 
body problem has a solution, regardless of the values of the 
free-streara Mach number and angle of attack. The analysis 
described in this Report can be easily started after solving the 
blunt body problem in a spherical, wind-oriented frame and apply- 
ing an auxiliary program to advance the solution through a spher- 
ical frame whose axis slowly rotates until it coincides with the 
body axis, t; simultaneously, the surface supporting the data at 
every step (which is a cone) slowly tends to become a plane, 
which will be the initial plane for the present computation. 

The second assumption is very convenient if the pointed 
nose has a circular cross-section and merges smoothly into the 
rest of the vehicle; in addition, the solution is accurate in the 
vicinity of the nose only if the angle of attack is zero. If the 
angle of attack is different from zero but less than the serai- 
aperture of the cone, a reasonable approximation can still be ob- 
tained in the vicinity of the apex. The namencal results, in- 
stead, can be grossly inaccurate if the angle of attack is larger 
than the serai-aperture of the cone. 

The assumption of a hollow intake in lieu of a pointed or 
blunted nose provides a speedy way of initiating the supersonic 
shock layer about the vehicle by surrounding the sharp rim of 
the intake with a locally two-dimensional attached shock; of 
course, an attached shock solution may not exist, and this limits 
the usability of the device, particularly for very low Mach 
numbers and high angles of attack. 

In the next two Sections, details will be given about the 
way the second" and third assumptions are implemented. 
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14. Pointed circular cone 


We assume that the trace of the shock on a plane, normal 
to the axis of the cone, is an ellipse if the angle of attack, a, 
is different from zero. Obviously, the trace of the shock is 
circular if a=0. We consider such a plane as tangent to the unit 
sphere, whose center is the apex of the cone, and also assume 
that, within the shock layer, 0 is practically equal to tan 0 
(here, 0 is the azyrauthal angle in a spherical frame, r, 0, (j)) . 

Therefore in Fig. i, where the plane is represented, 0 is equal 
to the distance between the point to which it belongs and the 
origin, 0, which is the trace of the axis of the cone. In the 
same figure, V is the trace of the vector, parallel to 
through the apex of the cone, E is the center of the ellipse 
representing the shock, 5^ and 6 ^ are the angles between shock 
and body in the leeward and windward side respectively, A and B 
are the semi-axes of the ellipse, e is the angle between shock 
and body on the x-axis, and b is the semi-angle of the cons. 
Note that 


G = A (1-C^/B^)^'^^ - b (163) 
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At zero angle of attack, the conical flow is governed by 
the equations: 


du/de = V 


(164) 


2 2 

dv/de = - u + (u+v cotan0)/(v /a -1) 


(165) 


Here, the three velocity components in the spherical frame are 
u,v, and w, respectively and a is the speed of sound. If o=0, 
the problem is solved by assuming a value for the shock angle, 
getting u, v and a behind the shock from the Rankine-HugOhiot 
conditions, and integrating (164) and (165) for decreasing values 
of 0 until V equals zero. Since such a condition should occur 
when 0 = b(that is, on the cone surface) and that hardly happens 
for the assumed value of 6^, the latter is corrected and the pro- 
cedure repeated until the condition is satisfied. To assure ac- 
curacy, I generally divide the shock layer into 300 intervals, 
which allows a simple method, such as Euler's, to be used. 


If the angle of attack is not zero, (165) must be re- 
placed by 


2 2 

v„ = - u + (u+v cotan 0 + w,/sin 0 )/(v /a -1 ) 

0 i{) 


(166) 


The same iteration technique is used to integrate (164) 
and (166), assuming that w is constant across the shock layer 
between M and N and between P and Q. In addition, we postulate 
that 


( e + 6_)/2 = 6 ( 167 ) 

2 O 

regardless of the Mach number and of the angle of attack. From 
(163) and (16?) A may be expressed as a function of B and C, 

26 -6„+b 

A = 

„ 2 .„ 2 . 1/2 

(1-C /B ) (158) 

Consequently, only tv;o parameters, B and C, are left to define 
the ellipse. They are related to 6^ and 6 ^ through the equa- 
tions , 
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B = b + (S^+S^)/^ 
C = (5^ - d2>/2 


( 169 ^ 


To find approximate initial conditions for the conical 
shock layer at an angle of attack, we call s the value of e at 
the shock. From the equation of the ellipse, it follows that 


? 2 2 2 2 1/2 
s _ Cc oS(|) + B[(B "C )/A sin ({i+cos <l>] 

(B/A)^sin^()) + cos^it) 


(170) 


and 



(171) 


at the symmetry plane. From the Rankine-Hugoniot equations, 

w* = us^..^sin s - V, (1721 

9 <}'4’ 9 

where u is the component of the velocity behind the shock along 

the normal to the shock wave, pointing inwards, and is the 

9 

derivative with respect to 9 of the component of the velocity 
along the tangent to the shock wave in the plane of Fig. 3. The 
component of the velocity in front of the shock along the normal 
is 


V sin(s±a) (173) 

00 

where the upper sign holds for point Q and the lower for point N. 
From (173), e is obtained through the Rankine-Hugoniot equations. 
From the same equations, 

V. = V Cs^.cos a±(s.^cotan s - 1) sin a] (174) 

ffl <X> m mm 

We have now the equations which are necessary and sufficient to 
determine B and C. From a computational point of view, an itera- 
tion procedure is convenient. The value of C is initially as- 
sumed as zero and a value of B is guessed; A is computed from 
(168) and the corresponding values of w^ on the windward and lee- 

” (J) 

ward side are evaluated. Then the equations (164) and (166) are 
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integrated between N and M and between Q and P. Two new values 
of 6^ and 6 ^ result, from which two new values of B and C are ob- 
tained, and the procedure is repeated uhtil convergence is 
achieved. 


Once the shock shape has been determined, the Rankine- 
Hugoniot conditions provide all the physical parameters behind 
the shock at any point. The pressure on the body is linearly in- 
terpolated between windward and leeward side; a equals the 
tangent to the cone angle; S is taken equal to its windward value 
everywhere except at the leeward point. The cross-flow velocity 
on the body is defined by a parabola: 

2 ^ 

w = b(j) + Ciji + dij) (175) 


where 

2 

From P and S, the temperature is determined and then q . Final- 

ly. 


... 2 2 . 1/2 

n = wd/(q -w ) 


(176) 


The values of P, S, a and n are linearly interpolated between 
shock and body on each meridional plane except the windward and 
leeward plane, where they are obtaihed from the converged solu- 
tion above. Such a linear interpolation between shock and body 
is rather simplistic and it is a far cry froim the real solution 
to the conical flow past a circular cone at an angle of attack. 
More parametric studies of the latter could provide a better set 
of initial conditions. The current procedure surely shows its 
age [8]. 


15. Attached shock around a hollow intake 


If we assume that the body begins with a hollow intake of 
a prescribed cross section, whose lip produces an attached shock, 
the latter and the physical parameters behind it can be evaluated 
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as follows. 

Let fl be the unit vector, normal to the lip, in the first 
t=constant plane, fl the unit vector, normal to the body surface 
in the same plane, and n the unit vector, normal to the shock 
produced by the lip, again in the same plane. Note that fl is 
contained in the t=constant plane, but N and n are not. 

For fl, (iO) and (i1) hold. Therefore, 

= 1/v^ , = -(by/b)H^ (177) 


with 


= [l+(by/b)^]''^^ (U8) 

Since the shock is attached to the lip, the latter is initially 
the cross-seciton of the shock itself; therefore, the unit vector 
normal to the shock in the t=constant plane coincides with fl. Ne 
may write: 


n = n fl + n,R (ir9) 

o 3 

and note that 

+ n^ = 1 (180) 

o 3 


Since 


n. = n^H, , n^ = n^H_ (181) 

1 O 1 2 O 2 

n^ and n^ are the only two quantities to be evaluated, in order 
to determine the initial direction of n. 

The free stream velocity is 

^ = ui + vj + wlc=un + v (182) 

00 00 00 00 00 CO 

Note that 

u = V S sina , v = V 2 cosa , w = V cosa (183) 

CO 00 00 00 00 00 
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according to (lOi) and (104). Ooviously, 


u = 7^ .n = u n.+v n„+w n.. = n^(u H,+v H-)+w n, 
=n + w n. 

O 00 j 


and 


V = \T -u n = (u -u n H. )i+(v -u n H„)j+(w -u n,)ft 

09 0909 090901 OOOOQi^ OOOO^ 

The velocity behind the shock is: 

rt A A ^ OVA ^ 

V = ui+vj+wk = un+v 


with 


~ Y-1 ~ 2y 1 

u = ru +— V ~ 

Y+1 <» Y+1 h 


so that 

Au = u-u = — ^(y/u -u ) 

09 Y+ I ' 09 09 

The velocity behind the shock must be parallel to the body: 


= 0 


From (189) t (18D, (185) and (188) we obtain the condition: 
[u +Au.n H.]N. + [v +Au.n H„]N, + [w +Au.n,JNj = 0 

00 oil 00 0 2 2. 00 jj 

The equation above, by virtue of (180), contains the only 
nown, n^. By introducing the notations: 


E=AN3-C+wJ„j^J^, F=BN^-2wJ„j^N3, G=BN3-.C 


(190) can be written in the form; 


( 184 ) 


(185) 


(186) 


(187) 


(188) 


(189) 


(190) 

unk- 
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(l-n^)(D+Fn2)^=(£+Gn^)^ 


(191) 


or 


6 K 'I 2 . 
n^+bn^+cn^+d = 0 


(192) 


where 


b = A[2(DF+£G)-F^], c = A[D(D-2F)+£^], d 


A = 1/(F^+G^) 


= -AD 


(193) 


Note that (192) has the same form as (150a) in [10] and a similar 

2 

physical meaning. If there is only one real root for n^, no at- 
tached shock is possible (the root corresponds to a value on the 

outgoing branch of the shock polar. Fig. 4). If there are three 
2 -- 

real roots for n^, the intermediate one corresponds to the physi- 



cal condition of a weak, attached shock (point B in Fig. 4). The 

2 

computational program chooses the intermediate value of n^ and 
then computes n from (180), u froiTi (184), Au from (188); it 

O OO 

follows that 

= u^/y (194) 

and, behind the shock. 


- 41 - 



Attached shock around a hollow intake 


P = in (Yt^„ - 

S = P + Y In (u/u ) (195) 

00 

u=Au n,+u , v=Au n»+v , w=Au n.,+w 

I 00 c. ^ i ^ 

From (33) and (3^), written with n in lieu of N, wa see 

that 

Cj = c[f^-(Cy/c)f2]-Gd = Cc(n^f ^+n2f2)-Gn^]/n^ (196) 

A very small step, At, may be taken, and a new shock shape as- 
sumed in the 5-plane, given by 

c = b + c^At (197) 

where c^ is defined by (196). Tnen, values of P, S, w, a = u/w, 
n = v/w as obtained in (195) are uniformly distributed between 

shock and body at the new station, for every value of 7. To 

speed up the subsequent calculation, only four points for every 
value of Y are considered, one on the body, one on the shock and 
two in between. A slight correction may be needed, in order to 
make a at the body compatible with the new body geometry at t+At. 
We leave n unchanged and, after computing new values of 
N^, N2 and according to (30), define 

a = - (N^+nN2)/N^ (198) 

according to (189). The initial values of v at the body are then 
determined by using (100). 
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16 . Preliminary testa . Circular cones at no incidence 


In what follows, detailed descriptions of the geometry 
subroutines for each body geometry are given. For each body 
geometry, the computations performed are reported with pertinent 
comments . 

A preliminary set of tests was performed on circular 
cones at no incidence. Any cross-section of a circular cone is 
defined by 

+ y"" = s'" (199) 

with 

a = At , a^. = A , a^^ =0 (200) 

Circular cones are used primarily to test the logic and the cod- 
ing of the mapping subroutines, and to get an idea of the accura- 
cy and reliability of the integration technique. 

The latter is the object of a first series of tests. The 
free-stream Mach number is set equal to 6, the angle of attack 
equal to 0 and the cone angle equal to 30^ (A = .57735). In Run 
1, the pointed nose model is used, and the computation starts at 
t=1 with 12 radial intervals and 2 circumferential intervals (of 
course, at no incidence, all meridional planes should contain 
identical results and all 7-derivatives should vanish) . No 
stretching is used. One mapping function is used, with h^^=a/2 
and 6=1. The resulting mapping is an identity. Runs 2, 3, and 4 
are similar, but the X-coordinate is stretched, with 
0=0.5, 1 and 2 respectively. In all four cases, the calculation 
is carried on for a large number of steps (100 for the first 
three cases and 200 for the fourth) and the differences between 
the final and the initial values of pressure (multiplied by 
10,000) are reported in Fig. 5. As the figure shows, the stabil- 
ity and accuracy of the computational technique is indeed remark- 
able (at least, when only one space-like dimension is involved). 



Preliminary tests. Circular cones at no incidence 


The errors are almost in round-off error order range. It must be 
noted that Fig. 5 is intended to give an appreciation for the 
order of magnitude of the errors; from step to step, they oscil- 
late in a wavelike form between body and shock, but they always 



remain confined to the ranges shown in the figure. 

These results also show that one can use rather strong 
stretching devices with confidence. Note, in Fig. 6 , the distri- 
bution of nodal points between body and shock in the four cases, 
and the very strong accumulation of nodes near the body resulting 
from the use of a stretching parameter equal to 2 . 

A second series of runs was made to test the mapping 
subroutines, using the same physical and geometrical inputs as in 
the first series. In the first run, two mapping functions are 
used in a sequence, the first with h^^=a /2 and the second with 
h 2 ^=( 1 +i)h ^1 . For both mappings, 6 ^= 62 =!; therefore, both map- 
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a = 2.0 

11-1— I I 1 I 


. .I-Q . . 

« t » » 


-i t_ 


,0 


-J 1_ 


Fig. 6 


pings are identities. 

In the second run, the effect of non-identical mappings 
is tested. Two mapping functions are used, with their hinge- 
points on the same location (h^^=h2^=a/2) . The values of 6 are: 
6^=l/2, Therefore, the first mapping transforms the circle 
into an ellipse, and the second transforms the ellipse back again 
into the original circle. 

In the third run, the effect of non-identical mappings 
and of variable values of 6 is tested. Here, 
6^ = .2+6(t-.2), 62=1/6.,. 

A fourth run was made to test whether the results were 
affected by making use of the symmetry of the flow with respect 
to the x-axis and, conse^quently, computing only in one quadrant 
of the z-plane. 

In all of the runs above the results were identical to 
the ones obtained in the first run, to within 4 significant di- 
gits. 


17. Circular cones at an angle of attack 


Next, flows past circular cones at an angle of attack 
were computed. For all of them, one mapping function is used as 
in Run 1, with h.,.,=a/2 and 6=1., The first case considered has a • 
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Circular cones at an angle of attack 


10^ cone in a free stream Mach number of M=2 at an angle of at- 
tack of 5°. The cone starts at t = 1 with a hollow intake; the 
initial mesh has 3 radial and 12 circumferential intervals; the 
number of radial intervals is doubled at t=1.5 and again at t=2. 
The stretching parameter is equal to 1 . The computation is con- 
tinued for 500 steps, bill t = 1880. At that step, the maximum 

value of c. . is less than 10“ and we can consider that the flow 
tt 

is entirely conical, A sequence of isobars plots for this compu- 
tation (Run 11) is shown in Fig. 7. 

Another calculation for the same geometry, Mach number 
and angle of attack was performed using a stretching parameter 
equal to 1.5; the radial intervals, initially in number of 3, 
were doubled at t=1.5 and t=2; the circumferential intervals, in- 
itially in number of 12, were doubled at t=2.5. A sequence of 
isobar plots for this case (Run 152) is shown in Fig. 8. At step 
500 (t=3^.4), c^^ does not exceed 0.0002. In Fig. 9 the circum- 
ferential pressure distribution on the cone at the final step of 
both runs is compared with the pressure distribution predicted by 
Jones C9]. In Fig. 10 a similar comparison is made for the radi- 
al pressure distribution at the windward and leeward side, as 
well as at 9=90*^. The agreement is very good. 

The evolution shown in Figs. 7 and 8 deserves a careful 
examination, since it shows how the flow, which is initially lo- 
cally two-dimensional and practically independent on the radius 
(as it should be in the proximity of the intake lip) gets twisted 
into a complicated conical pattern. In the first phase of the 
evolution, a pressure 'valley' and a pressure 'hill' appear; the 
former near the shock and the latter near the body (Fig. 11); 
further ahead, the valley disappears and the hill gets stabilized 
into the asymptotic pattern. 

Another shockless flow field is obtained, for M=2, about 
a 5^ cone at 10*^ angle of attack (Run 121). The computational 
parameters adopted for this case are the same as for Run 152; 
only the stretching parameter has been set equal to 2. A se- 
quence of isobar plots is shown in Fig. 12. The asymptotic pres- 
sure coefficient, as computed, is compared in Fig. 13 with values 
obtained by other Authors [11]. Note the formation of a pressure 
valley on the leeward side of the cone, anticipating the forma- 
tion of a cross-flow shock, which occurs for higher incidences. 
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18 . Elliptic cones 


Any cross-section of an elliptic cone is defined by 



1 
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Elliptic cones 



with 


a = At , = k 

b = Ba , = BA 


( 202 ) 


The derivatives of F, as defined by (35), are 

F = 2x/a^ , F = 2y/b^ , F = - 2A/a (203) 

X ’ y ’ t 

To obtain a perfect circle as the image of the body in the 
C-plane, one mapping function can be used, with 

f, q2,1/2 * _ o2J/2 . _ 

li,|^ — sCl— B ) f ~ A\1“B ) I o ,j — 1/2 

The other two points of interest are the lower and upper inter- 
section of the body with the x = 0 line, viz. 


= - ib , 1^31 = ib 


- 50 - 



Elliptic cones 


•4h 



from which it follows that 

• • 

hg^ = - iAB , h^^ = iAB 

A computation was performed (Run 16) with A = 0.6, B = 

0.5, = 6 at an angle of attack equal to 16°. The initial lip 

is assumed to be at t=0.2. The computational grid begins with 3 
intervals in the radial direction and 12 intervals circumferen- 
tially; the radial intervals are doubled at t=0.ij and again at 
t=0.6. The circumferential intervals are doubled at t=0.6 only. 
The stretching parameter, a, is made equal to 1 . In Figs. 14 and 
15 isobar plots at selected steps are shown. After 1500 computa- 
tional steps, the station defined by t=12.37 is reached and the 
computation is halted. At this stage, the computational grid ap- 
pears as in Fig. 16. In Figs. 17 and 18 curves of constant P, S, 
o and n are shown. In all these figures we may note an indenta- 
tion in the outer shock at the leeward symmetry point. Does it 
indicate a programming error? A careful analysis of outputs of 
preceding steps led us to the conclusion that slight inaccuracies 
may result because of insufficient resolution, but that a similar 
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Fig. 11 


trend is justifiable on physical grounds; the indentation is ex- 
pected to disappear after a large number of additional computa- 
tional steps. In other words, the flow is not yet conical on the 
leeward side after 1500 steps. Let us examine then how the shock 
shape varies from t=0.2 to t=12.37 and how the conical flow pat- 
tern sets in. 

At the beginning, we assume that, at every point on the 
lip, a shock is formed according to the analysis of Section 14. 
Such a shock is evidently two-dimensional at each point since no 
signal from neighboring points on the lip can be transmitted la- 
terally due to the supersonic nature of the flow. Therefore, the 
initial shock and the consequent distribution of values between 
shock and body are not conical. The rate of change of the shock 
in the physical space is controlled by the values of c in the 
mapped plane for every cross-section. We can, thus, analyze the 
evolution of the distribution of c as a function of 9, as t in- 
creases. In a conical flow, c^^ should be constant in t for every 
point. In Fig. 19, computed values of c^^ are shown at selected 
steps. Initially, such values are very-high near 0=0 where the 
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Fig. 13 


curvature of the cross-section of the body is the highest and, 
therefore, the deviation of the initial two-dimensional shock 
from the conical one is also the greatest. The windward side of 
the flow, however, tends to get stabilized very rapidly, thanks 
to rather low Mach numbers in the shock layer (of the order of 
2). All perturbations tend to be carried towards the leeward 
side since the velocity vectors are so oriented (n>0). The shock 
tends to move outwards for values of 0 greater than 0, and the 
point of maximum c shifts rapidly towards the leeward symmetry 
plane. Note, incidentally, the refinement in the curve of c^ 
which occurs after the number of grid points in the circumferen- 
tial direction is doubled. 

As t increases, we expect the peak in the c^-distribution 
to move further towards 9=iv/2, and eventually reach it. At this 
stage, however, the locus of the peak is shifted only very slow- 
ly; for an observer watching the cone from above, the angle under 
which the peak approaches the leeward symmetry plane is about 
10®; that is, the numerical perturbation travels along a charac- 
teristic. A finer mesh could show a better behavior of the shock 
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shape in the vicinity of the leeward symmetry plane, but we do 
not expect a perfectly convex shape to appear, even after a 


I 


Elliptic cones 



RUH» 16 H» 6.60 K=1501 T« 

. DREF«0.8290 HAX«0.66e8 


4 


12.3728 




I 

I 


I 

\ 



fiUHs 16 Ms 6.68 K>1581 T« 12.3728 


0REF>8.8268 KAX>e.3868 5 


Fig. 18 


prohibitive number of additional steps. 

Quite reliable results on the same case have been ob- 
tained by Pandolfi [171. We present some comparisons of our 
results at the last computed step with the latter's in Figs. 20 
through 22. The shock shape is compared in Fig. 20 and the pres- 
sure distribution in Figs. 21 and 22; Fig. 23 shows the 
n-distribution along the body. We note some rather large 
discrepancies in the upper portion of the body or, more precisely, 
from the conical 'stagnation point' upwards; the stagnation point 
itself, however, coincides in both calculations. The value of n 
at the body is, in our case, obtained from (101), which in turn 
depends on (100), (99) and (98), The latter is a sensitive equa- 
tion and small changes in its coding may produce large variations 
in results, for flows such as the present leeward side flow. 
Pandolfi and myself have commented on this difficulty in a previ- 
ous paper [5], but this does not imply that his elliptic cone 
program and the present, general-purpose program are similarly 
coded! In particular, I believe that Pandolfi did not compute v^ 
as I do and this is sufficient reason to account for the 
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Fig. 19 


discrepancy. At this time, it is hardly possible to make a 
statement about the accuracy of either result. 

A second computation was performed (Run 24) with A=.3640, 
B=.1667, M^= 2, at an angle of attack equal to 10®. The hollow 
intake is assumed at t=1; 24 intervals are used circumferential- 
ly; radially, we start with 3 intervals, doubling them at t=1.25 
and again at t=1.6; the stretching parameter is taken equal to 1. 

The computation is halted at step 1000 (t=11.112). A se- 
quence of isobar plots is shown in Fig. 24. In Fig. 25, the 
values of P around the body are compared with results by other 
Authors [12,13]. The agreement between the present results and 
those obtained by relaxation is excellent in the windward side. 
On the leeward side, our results seem to overshoot the expansion 
and the 'shock' appears too far to the left. 

In order to assess the reliability of the calculation and 
the effect of resolution in the radial direction, the calculation 
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Fig. 23 
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has been repeated increasing the stretching parameter to 1.5 (Run 
124). To avoid difficulties in the initialization, the initial 
conditions are taken from the solution of a circular pointed cone 
at zero angle of attack, the cone semi-angle being equal to 20^. 
The initial cross section is located at t=0.1; doubling of the 
radial intervals occurs at t=0.125 and t=0.l6. In Fig. 26, our 
results after 1500 computational steps are compared again with 
those of [ 13 ]. The overshoot is slightly improved, but the mesh 
is definitively too coarse circumferentially to resolve the shock 
properly. Whether the shock could be 'captured' in the right po- 
sition by using a finer mesh, still remains an open question. 
One thing which appears from the evolution of the calculation is 
that the improper collocation of the shock is closely related to 
a progressive deterioration of the flow downstream of it. Some 
wiggles develop in the pressure distribution on the leeward side 
of the body, as t increases. In conclusion, fitting of the im- 
bedded shock is necessary for accuracy on the leeward side, at 
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least with the present resolution. 

Finally, the flow past an elliptic cone with a 1:14 axis 
ratio was evaluated (Run 25). in this case, A=0.3640, 8=0.0714, 
M^=2, at an angle of attack equal to 5°. The computation is 
started at t=0.1, assuming, as in Run 124, that the inital condi- 
tions are given by the solution of the conical flow at zero in- 
cidence about a circular cone having a diameter equal to the ma- 
jor axis of the elliptic cone. The number of intervals and their 
doublings are as in Run 124. The computational grid and a set of 
isobars at step 2000 are shown in Fig. 27. In Fig. 28, the pres- 
sure coefficient is compared with the results of other Authors 
[12,133. In this case the agreement is better, despite the ex- 
treme coarseness of the mesh. This case is, in my opinion, a 
good example of the efficiency of the computational technique. 
Note that Marconi's results have been obtained with a method very 
similar to the present one, but not employing the X-scheme, and 
with a careful fitting of the imbedded shock; note how much finer 
his mesh is. Nevertheless, our present 'shock' is in the compu- 
tational cell where the fitted one appears; about the wiggles and 
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T« 



Fig. 27 


general discrepancy of results behind the shock, we can repeat 
the comments to Runs 24 and 124; we may also note explicitly that 
our 'shock* is isentropic. We are sure that no higher resolution 
is needed for the calculation of these flows, once the code is 
supplemented by a shock-fitting routine. 




1.8787 
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19 . Cambered wings 


To test the program on shapes which are of current in- 
terest, without making the mapping unnecessarily cumbersome, we 
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29 


generate a cambered wing as follows (Fig. 29). In the t-plane, a 
circle centered at the origin and with a radius D is considered 
as the image of the camber line of the wing. Another circle, 
concentric to the first, with a radius 1, is the image of the 
wing itself. The -plane is obtained from the C-plane by a 
simple translation: 
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Z2=t+io , i=z^-ia C204) 

Let be the point of the c-plane whose image in the Z 2 ~plane is 
h. Let the Z-plane be obtained from the Z 2 ~plane via a Joukowski 
mapping: 


Z = 



(205) 


In the Z-plane, the camber line becomes a circular arc counted 
twice. Finally, let a second translation be applied to get the 
z-plane from the Z-plane: 


z=Z-i3 , Z=z+iB (206) 

We want the cambered wing to have a maximum thickness equal to 2, 
with its points z^ and Zg equal to -i and +i, respectively. Let 

C = B-(B^-4h^)^^^ (207) 

where B is a prescribed constant. Clearly, C is the distance 
between point P in the Z-plane and the origin, if B is the ra- 
dius of the camber line. Therefore, 

.2 

“ ^ «08) 

Let us prescribe the variation of x^ with t, for example by say- 
ing that 


= At (209) 

where A is a prescribed constant. It follows that 

h = At/2 (210) 

Since, as it can be seen in the Z 2 -plane, 

= h^+o^ (211) 

(208) and (211) provide D and a as functions of h: 
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D = (H+h2/H)/2 , 0 = (M-h^/M)/2 (212) 


where 


M = I [C+2BO^^^] (213) 

and, because of (210), of t. The remaining two unknown parame- 
ters, 1 and B, are determined by the conditions imposed upon 
points Q and S in the z-plane: 


1 ^ =-l " ^ 

" ^ ^ P 

Therefore, 1 can be determined by solving the third degree 
tion (which has only one real root): 

(1-1) (o^- 1^^ +h^l = 0 

and B follows from (215). The t-derivatives are: 


(214) 

(215) 
equa- 


(216) 


\ 

DMC^. + hh^(C-2o) 
°t " (2M - C) M 

D^. = (hh^ + oa^)/D 
DID -0O. 

1=2 

(31-2) 1 - D^ 


(217) 

(218) 

(219) 

( 220 ) 



(at . V (1 




_hi_. 

( 0 . 1)2 - 0.1 


( 221 ) 


Since we do not have an implicit form, F(x,y,t)=0 for the 
wing contour, and we need F^ and F^, we may proceed as fol- 
lows. First, we use (10), (15) and (22) to evaluate 
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Cambered wings 


V“x^ ^ + 2(^2^ (222) 

2y = -2(0y-(Oy) = - ^ [2(1-(j)^) + £^ 2 ^ + Y'^2 “ ^x ^ T 

(note that £ + iS, defined by (9)t is not an analytic function; 
therefore, £^^ S^) . Now, since the image of the wing in the 
C-plane is a perfect circle, b^ vanishes identically. Therefore, 
(38) tells us that 


F r KS , = K£ (224) 

where K is a function to be determined. By writing that F =F 

, » xy yx 

and using (21) we see that K is defined by 


- KyC = K ^ <^2 (225) 



Fig. 32 


The equation is satisfied by K = G, as (225) proves. Therefore, 
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“ . 

fy 

= G& 

(226) 

In addition , starting from (41) and 
and b = 1, we obtain 

the 

first of (31), with bY=0 

ft = - 

SFy) 

= 1^1 - "t 

(227) 

1 * 

As is evident from Fig. 

29 

and Eqs. (208) 

through 


(210), only one mapping function is needed to obtain the circle 
image of the wing. The pertinent value of 6 is obviously 1/2 and 
the position of the hinge-point in the physical plane is that of 
point M, that is, 


h^^ = At - iB , = A -iB^. (228) 

The other two points of interest are the lower and upper inter- 
section of the wing with the x = 0 line, viz. 

hg^ = - i , ^31 = i ^229) 

Their t-derivatives obviously vanish identically. 

Two examples are given. In the first (Run 28), the angle 
of attack is zero and the free stream Mach number is 3. The 
values of A and B are 0.4 and 5, respectively. The stretching 
parameter, a, is equal to 1. The hollow intake is at t=.5; ini- 
tially, we have 3 intervals radially and 12 circumferentially; 
the radial intervals are doubled at t=2 and 3.5; the circumferen- 
tial intervals are doubled at t=2. Some isobar plots are shown 
in Fig. 30; the computational grid and level lines of o, n and S 
at step 550 are shown in Fig. 31. Finally, in Fig. 32 the pres- 
sure distribution around the body at step 550 is shown. 

In the second example (Run 29). the angle of attack is 5^ 
and B IS 20. All the other values are unchanged. Again, we show 
isobar plots in Fig. 33, the computational grid and level lines 
of a, n and S at step 1100 in Fig. 34, and the pressure distribu- 
tion around the body at step 1100 in Fig. 35. 
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Cambered wings 



RUH* 23 «■ 3.00 K« 100 

OREF-0.1000 HAX>8.0000 2 , 

T» 4.2265 1 



RUH» 29 K« 3.00 K« 408 
CREF»0.10O0 MAX«0.9080 2 
T- 16.0556 



RUH« 29 «• 3.00 K" 550 
OREF-0.1000 nAX«1.0O00 2 
T» 22.5321 



RUN- 29 M- 3.00 K-U08 

'-OREF-0.1000 KftX*l 2000 2 

T» 35.3109 



RUN- 29 H- 3.00 K-1500 
DREF-0 1008 NAX«1.600a 2 
T« 40 5016 



RUN- 29 «- 3.00 K-2eO0 

OREF-8.1000 nAX-l.60Oe 2 , 
I- 42.9527 


Fig. 33 


- 76 - 



Butler's wing 



Another set of tests has been run in connection with work 
in progress at the University of Salford, UK. A body is defined 
as follows. For the first 20^ of its length the body is a circu- 
lar cone. The remainder of the body has elliptic sections which 
become more eccentric as the sharp trailing edge is approached. 
The values of a and b as defined by (201) are: 


a = At 


b = At 0<t<.2 

b = At [1 - •2<t<1 


(230) 


This body was considered in [14] and we will call it Butler's 
wing for brevity. From (230) it follows that 


a 


t 


A 
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Butler's wing 


b^. = A 0<t<.2 (231) 

‘'t = I - 5TC « - -2<t<1 

and 

= 2x/a^, F„ = 2y/b^, F. = -2(Ax^/a^ + b.y^/b^) (232) 

X y V w» 

The definition of the mapping and the location of the hinges are 
obviously similar to the ones given for the elliptic cone. 

The value of A was chosen as 0.29814; the free stream 
Mach number is 3.5. The first calculation (Run 14) was made at 
no incidence, with a=2, using 12 intervals radially and beginning 
with 6 intervals circumferentially; the number of circumferential 
intervals was doubled at t=.4 and again at t=.6. Typical isobar 
plots are shown in Fig. 35. In Figs. 37 and 38 some -experimental 
values [14,15] are compared with the results of our calculation. 
In Fig. 39 our results are compared with Butler's. It is obvious 
that our results are closer to the experimental values, which 
should not be a surprise since our calculations took advantage of 
twenty years of improvements in numerical techniques; the match 
is not perfect, however, and possible effects, for example of 
viscosity, remain to be explored. 

A similar comparison is given in Fig. 40 for the same 
wing, same free stream Hach number, at an angle of attack of 10°. 
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RIW* I< H= 3.50 K* 350 T« 0.6354 OREF'0.1000 RUh= 14 «» 3.56 <■ 400 T« 0.7219 O3EF«0.1030 

'•nftX»6 7C06 2 -- «flX»0.7O00 2. 

t ! 


i 

J 

i 


Fig. 36(a) 


i 
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Fig. 36(b) 
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Butler's wing 


□ REF. I 
O REF. I 


CALCULATED 


cn ^ 



Butler's wing 


EXPERIMENTAL 
REF. 15 


O 

A 

□ 

O 


t = 0.3 
t = 0.5 
t = 0 7 
t = 0.85 




Butler's wing 


PRESENT THEORY 
BUTLER I REE 14) 







21 . A simple fuselage-arrow-wing combination 


A configuration which is closer to a realistic airplane 
geometry with a distinctive differentiation between a fuselage 
and an arrow wing has been defined as explained in this Section. 
To maintain the analysis as clean as possible, the entire 
geometry has been defined analytically. The configuration has 
been chosen trying to minimize causes of formation of imbedded 
shocks . 

The top view of the airplane (Fig. 4-1 ) shows a straight 
'leading edge*, 


X = a = At (233) 

where A is an input value. The nose of the airplane is an ellip- 
tic cone, whose cross-sections are defined by 



with 


„ _ t(A+j5t) 

® ~ ~TW~ 


where 


t = k - 2/t, 


^ = - M< 


U J ^ - I / 

and is an input value. It follows that 


B = 1 + t(}$-A^) 
^ ^ (l+?!t)^ 


(234) 


(235) 


(236) 


(237) 
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A simple fuselage-arrow-wing combination 


and B = 0, = A at t = 0 and B = 1 , B^ = 0 at t = t^. 

Another line, 


X = Xg(t) (238) 

defines the 'center* of a wing which begins at t = t^, the latter 
being an input value. The wing is a Joukowski profile, modified 
by an additional thickness to give it some structural strength. 
The basic profile is obtained by mapping a circle, 

T = g + 1 (239) 


onto the z-plane, via the transformations 

+ y^ = + T + k^/F , z = x^ + iBy^ (240) 

The additional thickness is obtained by adding a value, y, given 
by 


= e(x-: 


TE 


) 

^TE"^ 


(241) 


where x^g is the abscissa of the trailing edge of the Joukowski 
profile, to each y obtained from (239) and (240); B is a function 
of t, to be defined later. Beyond t^, an elliptic fuselage, 



(242) 


begins to be differentiated from the wing, until at t = t^ (an 
input value) wing and fuselage become detached. The end of the 
airplane is at t = tg (another input value). 

To assure a smooth transition from the nose to the wing- 
fuselage arrangement, the Joukowski profile must be made to coin- 
cide with the elliptic cross-section at t = t , and no additional 

o 

thickness must be used there. This is obtained by imposing the 
conditions: 
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A simple fuselage-arrow-wing combination 


g=e^=0, l=l^=(a^+l)/2, k=k^=(a^-l)^/^/2, 0=B^=O (243) 

at t = t . 
o 

Since the leading edge of the profile corresponds to c=g+l 
and it must be at x = a, we have the condition: 


a = x + g+ l+ —rir 
0 ^ g+1 


The trailing edge, x = x^^, corresponds to c=g-l: 


(244) 


= ^0 e - 1 + (245) 

At t = t^ , we want 

(x^.^)^ = 1 (246) 

and we also want to have the cusp of the profile at x=1 : 


- 1, - g, 


(247) 


Finally, we want 




( 248 ) 


where is a prescribed value, at t = t^ and constant for all 
'*'2.^1* Note that, at t = t^ (245) can be replaced by a combina- 
tion of (244), (245), (246) and (247): 


41 


1 l^+g. 


(249) 


We should also prescribe g^ , in addition to A, t^, t^ , t2» t^, 
and B^ ; the two parameters, g^ and B^ control the thickness of 
the wing. Then, (249), (247) and (244) yield 1^, k^ and » 
respectively. We will also impose that 


Xq = E (t - t^) 


( 250 ) 
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A simple fuselage-arrow-wing combination 


between t^ and t^. Therefore, =■ E(t^-t^) and it is clear 
that E and t^ cannot be prescribed independently. We choose to 
prescribe E. 

Between t^ and t^ , g, Q and 1 can be interpolated linearly, 
and is given by ( 250 ). Therefore, k and are obtained from 
(244) and (245). In particular, 

k = [(a - x^ - g - l)(g + 1)]^/^ ( 251 ) 

Between t and t., one must also determine the point of intersec- 
tion between wing and fuselage. The wing is defined by ( 239 ) » 
(240) and (241); letting X = x-x^, T = g + 1 cos (t> , and noting 
that 


(g + 1 cos 41)^ + l^sin^(j) = 2gr + 1^- g^ (252) 

the Cartesian coordinates of a cross-section of the wing are 
given by 


X = r (1 + 


2gr+i^-g^ 


(253) 


y = 1 sin 41 (1 - 
From the first of ( 253 ) 


2gr+r-g‘^ ""te ® 


^ -d^-g^^k^-2gX) ^ [(l^-g^+k^-2gX)^+8gX(l^-g^^^ 

4 g 


1/2 


cos 4 > 


r-g 


r-gN2i1/2 


sin = [1 - (-^^) ] 


(254) 

(255) 


and y can be obtained from the second of ( 253 ). The values of x 
and y which satisfy the latter and ( 242 ) simultaneously, which we 
will call Xjj^, y^jj, define the wing- fuselage intersection. 

Beyond t^ , ( 250 ) is no longer valid. The trailing edge 
abscissa is defined instead: 
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A simple fuselage-arrow-wing combination 


2 '*^■'*^1 

■ ' * (*‘ 2 -’ ) 7 "='==^"tpt7 


( 256 ) 


and g: 


. t+to-2t^ 

(a-x^g) sin (| — ) 


(257) 


Here, (244), (245), (247) and (248) are valid; therefore, (244) 
can be replaced by a relation similar to (249): 


® ■ ^TE 1+g 

which yields 1; then, k is obtained from 

k = 1 - g 


( 258 ) 


(259) 


and, from (245), 


^0 " ^TE ^ ^ 

A view from the top of a typical geometry is shown in Fig. 
41. Some typical cross-sections of the airplane are shown in 
Fig. 42. They have been obtained with A=0.5, B=0.05, E=0.3195, 
g=0.035, , t^=21.3287, and t2=31.9930. and g^=0.6. 
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22 . Mappings for the arroi^-wing airplane 


In the first section of the airplane (t < t^), where the 
cross-section is an ellipse, only one mapping function is needed, 
and it is defined as for the elliptic cone. 

In the second section, we need three successive mappings, 
one to open up the wing and two more to eliminate the corners at 
the wing-fuselage intersection. The first mapping is similar to 
the one used in the first section. The radius of curvature of 
the leading edge of the wing is 


l[2(g+l)-(a-Xp)]^ 

U-Xp)(g+1) 


The first hinge is defined hy 

h^^ = [a(a-p)]^/^ = 1/2 

The second and third mapping are defined by 


^21 ^IN ~ ^ yiN» ^2 ®IN 

^31 " ^IN ^ ^IN’ *3 " ^IN 


The intersection points, are determined using a trial- 

and-error routine. Eqs. (254) and (255) are used. In the code, 
the symbols; 

= 1^ - g, d = - 2gX, = 4gr + d 

A = 2gr + D^, E, = 1 - k^/A, 

~ ^^IN ~ ~ a), Ej = 1 sin<^ , E^ ~ ~ ^TE^ 

are used for simplicity. Then the matching of 
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Mappings for the arrow-wing airplane 


y'i’ - E, Ej + eI 
(which is the second of (253)) and 

which is ( 242 ), defines and The t-derivatives of Xjj^ 

and are obtained by differentiating the quantities above: 

■’t ■ “it * ("a^t - * ='t)s”2t - * ■■*) * 

r* = [-rdj ♦ D,(xg^/g + x^) * XDj^l/Dj - rg^/g 

(cos4)^ = (g^ - - l^cos4')/l, (sin<t>)^ = - cos4' (cos<l> )^/sin4' 

^t " 2(rg^/g + r^)g + 

E,t - [-(k^)t + (1 -E,)Aj]/A 

= l^sin<!> + l(sin4')^ 

^4t " ^t^^IN “ ^TE^ ~ ^^TE^t^ 

" ®1®3t '*' ®3®1t ^"^^4t^2 '*' ^®4®2t^®2 

Obviously, (yjjj)^ and Xjj^)^ must be obtained by a second trial- 
and-error procedure, the object of which is to match the value of 
(yijj)t above with the derivative of (242); 

Once all these values are found, the value of 5 2 and 6^ is found 
as follows: 

^ ^/^2» (sin<t>)^ = r^cos't'/E^, = 2gf^ 

®1x ■ ('-’^1 V''- E2^ - l/(xjj.-a) 

- 92 - 



Mappings for the arrow-wing airplane 


Ej = Ej^E, * EjE,^, Eg - E^^E^ * 2E^E2^_ T, - Eg h- EjEg 

2 

Dj = arctanT^, = arctan(-B 

«2 = ^3 = 

To compute the derivative of 6 2 with respect to t, we use 

the following expressions: 


r 


xt 


(‘t- ^»2t 


)/»2 


(sin4>)^^ = [r^^cos(Ji + r^(cos(l))^- (sin()) )^E^^]/E^ 

E,^t = [0-E,)4,t- B,gA^- E,^SJ/A, = - e|J(^jj,)j- A] 

®3xt = ®4xt = ®t 


%t ®3xt^1 ®3x^1t ®3t^1x ■" ®3®1xt 


®6t “ ^4xt®2 ■*■ ®4x^2t ■" ^^^4t®2x ®4®2xt^ 


D 


®3t = ^t ^ ®2t®6 " ^2^6t)/^^ " " 3 ) 


*^2t ^3t ■ “ ®2^®4t " ®3t^/’^ 


(262) 


In the vicinity of t^ , and become larger than 2 and 
then tend to 2 as t tends to t^ ; in addition, their t-derivatives 
tend to infinity. Therefore, it is convenient to replace (261 ) 
and (262) by the following expressions: 


5^ = 2 +t(B 2 +BjT), 6^^ = B2+2 B^t (263) 

The symbols in (263) are defined as follows. Let 5. and te 

xo xxo 

the values of d. and 5,. at a cross-section where 6. becomes 
1 it 1 

greater than 1.85; let t* be the value of t at such a cross- 

section and 


L 
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Mappings for the arrow-wing airplane 


( 264 ) 

The above fit merges smoothly with (261 ) at t = t* and goes 
smoothly to 5^ = 2 at t = t^ . 


The ether two points of interest are the lower and upper in- 
tersection of the fuselage with the x=0 line, viz. 


h 


41 


= - ib. 



= ib 


23 . Additional formulas used beyond _t^ 


Beyond t^ , the following formulas are used to get the first 
and second derivatives of the geometric parameters: 


4-| - (t-t^ )/( t2~t^ ) , ^1t ~ ^/(t2~t^) 

d = arcsin A^, d^ = A^ ^/(1 -A^)^ 

A 2 = 2(At2-1 )/r 

A^ ~ ir A^ ^ ( t+t2~2t^ )/2 , ^55 ~ ^^1t^^ 

^TE ^ ^■*■'^^2’ ^^TE^t ^2^t 
d^ = (a-Xjg)/2, = [A-(x,pg)^]/2 

g = rd^sin A^, = d^^g/d^ + Td^A^^cos A^ 

‘4 ■ ^ 4 t ■ [d,<i,t*'>(gait*'>iet ’]/‘4 

l-(d,t6^)/4, Ij - (d,t+64^)/4 
k . 1 - e, k,. . 

^0 ^ ^ot ^ ^^TE^t"^^^t 

X = (l+g)k, = 2(ll^-gg^) 

For the second and third mappings, 


(265) 

(266) 

(267) 

(268) 

(269) 

( 270 ) 

(271) 

(272) 

(273) 

(274) 

(275) 

(276) 
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24 . Some results on fuselage- arro w-wi ng combinations 

Ve will report here the results of seven runs made with the 
geometry above. Runs 50,51,61,62,63 and 64 were made with the 
following values of geometrical parameters: A=0.5, B=0.05, 
E=0.3195, g=0.035, t^ =21. 3287, and t2=31.9950. Run 52 was 
made using the following geometrical parameters: A=0.25, B=0.05, 
E=0.1831, g=0.06, t^=2n, t^ =27.81 93, and t2=38.3906. The angle 
of attack is equal to zero in Runs 50, 51 and 52, equal to 3 de- 
grees in Runs 61,63 and 64, and to 6 degrees in Run 62, The Mach 
number equals 2.36 in Runs 50 and 64, 2.96 in Run 63, and 4.63 in 
Runs 51 ,52,61 ,62. 


In the following Pigs. 43 through 62, computed values of C 

p 

are compared Vfith measurements made at NASA (Ref. 16). The 
values of X quoted in the figures are, in the notation of the 
present Report, values of t/t 2 * Similarly, the values of Y are 
-values of 
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Fig* 45 
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RUN 61 
Y*0.0467 
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Fig. 53 

; 
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Fig. 54 
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RUN 63, 

■ Y=0.1667 



Fig. 57 


ami 
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Fig. 59 
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Fig. 60 
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